Cluster algorithms for general-S quantum spin systems 
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We present a general strategy to extend quantum cluster algorithms for S — j spin systems, 
such as the loop algorithm, to systems with arbitrary size of spins. In general, the partition function 
of a high-S* spin system is represented in terms of the path integral of a S = | model with special 
boundary conditions in the imaginary-time direction. We introduce additional graphs to be assigned 
to the boundary part and give the labeling probability explicitly, which completes the algorithm 
together with an existing S = ^ cluster algorithm. As a demonstration of the algorithm, we simulate 
the integer-spin antiferromagnetic Heisenberg chains. The magnitude of the first excitation gap is 
estimated as to be 0.41048(6), 0.08917(4), and 0.01002(3) for S = I, 2, and 3, respectively 



PACS numbers; 75.10.Jm, 02.70.Lq, 05.10.-a 

O ' The world-line quantum Monte Carlo (QMC) method is one of the most powerful tools in numerical investigations 
^ ' of quantum spin systems Q. One of the main advantages of the QMC method over other numerical ones, such 
as exact diagonalization and density matrix renormalization group (DMRG) method, is that it is applicable to 
^ rather large systems in any dimensions and can estimate their physical quantities statistically exactly. However, the 
^ ■ conventional algorithm, based on local updates of spin configurations (world lines), suffers from strong correlations 
c/2 I between successive configurations at low temperatures or in the vicinity of a second-order phase transition. The 
4-^ ■ diverging auto-correlation time virtually makes simulations slower and slower, and finally it becomes practically 
Jh ' impossible to simulate larger systems at lower temperatures. This drawback is called critical slowing down. 
Ch Recently, the inventions of the loop algorithm |^,^ and of its continuous-time variant Q have led a great improvement 

of the QMC techniques for the S — ^ XXZ model 1^. The loop algorithm, which is a kind of cluster algorithms, 
(~| realizes updates of the configuration by flipping non-local objects, referred to as loops. It has been shown that it 
O is fully ergodic and drastically reduces the auto-correlation time, often by orders of magnitude, especially at low 
I temperatures. Furthermore, by using the continuous-time version of the algorithm, one can completely eliminate 

[ discretization error originating from the Suzuki- Trotter decomposition; simulations can be performed directly in the 
. so-called Trotter limit. 

In general, it is a highly non-trivial task to construct an efficient cluster algorithm for a given system, because 
symmetry or special properties of the target system should be taken into account explicitly in its construction. As for 
^ the spin systems, the development of cluster algorithms for higher-S* models remains as an important and challenging 
problem. A cluster algorithm for the general-S" XXZ model in the discrete-time formulation has been proposed by 
Kawashima and Gubernatis |6|. Unfortunately, their algorithm is rather complicated (105 different graphs appear 
even in the S — 1 case) , and moreover the Trotter limit is not well-defined in the algorithm. Another kind of algorithm 
in the discrete-time representation has been used in Ref. [Q. 
. More recently, Harada et al. proposed a continuous-time loop algorithm for the 5=1 antiferromagnetic Heisenberg 
model in which the S — \ system is mapped into a path integral of an 5 = ^ system with special boundary 
conditions in the imaginary-time direction. In the present letter, we generalize their method to construct cluster 
''O ' algorithms for systems with arbitrary size of spins. For simplicity, we consider a spin-S* antiferromagnetic Heisenberg 
C . model on a bipartite lattice as an example. Generalization to other models is straightforward as discussed later. 
We consider the Hamiltonian of spins and A'b bonds, defined as follows: 



^ = - L i^^^l + ^tS'l - St SI) , (1) 

where Sf {a = x, y, z and i = 1, ■ • ■, iVg) is the a-component of the spin-S* operator at site i. 

In order to construct a cluster algorithm for 5* > 5, it is crucial to represent each spin operator Sf as a sum of 
5=5 Pauli operators {cf^}. Following Ref. ||], we substitute Sf = 5 '^f ^^^'^ the Hamiltonian (|l|), which 

yields the following Hamiltonian of 2SNs spins and 4S'^A'b bonds: 

-, Wb 25 

^ = -i E E «.-^ + <A - • (2) 
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FIG. 1. Six graphs for the local boundary weig ht for 5=1. The upper (lower) three circles denote subspins at r = 
In general, (25)! graphs can be labeled uniquely in terms of corresponding permutations of 25* integers as shown in the bottom 
row. 



We refer to the S = ^ spins in Eq. (||) as subspins hereafter. In terms of the Hamiltonian (||) defined on the extended 
phase space, the partition function of the original Hamihonian (|^) can be expressed as 

Z = Tr(exp[-/37i:]P) . (3) 

Here, we introduce a projection operator P, which is a direct product of local symmetrization operators {Pi} {i — 
l,---,Ns). Each Pi acts on 2^'^-dimensional Hilbert space spanned by {crf^} {iJ, — 1,---,2S), and projects out 

unphysical states with Sf < S{S + 1). Note that P commutes with H, because H is invariant under the exchange of 
subspin indices at each site by definition. 

Applying a Suzuki- Trotter decomposition for the exponential operator in Eq. and inserting complete sets of 
eigenstates of {crf^} between the exponential factors, we obtain the following path-integral representation of the 
partition function: 

Z = Y,W{C)P{dC) , (4) 
c 

where C = {Ci^^^r} (* = 1, • • • , ^s, A* — ■ ■ , 2S', and < r < /3) denotes a world-line configuration of subspins 
and dC — {Ci^^^OiCj^A'^/j} (i — l,---,iVs and /i — l,---,25) denotes a configuration of subspins at the boundaries 
in the imaginary-time direction. We take the continuous-imaginary-time representation, and Ci^^^r denotes the /i-th 
subspin direction (-1-1 or - 1) at i-th site and imaginary time r. The weight P{dC) — Yii P-ii^^i) originates from the 
projection operator P, which can be interpreted as soft boundary conditions in the imaginary-time direction. Each 
(4S')-body local boundary weight Pi{dCi) takes a value of the inverse of the number of different configurations which 
are connected with dCi by permutation operations, i.e., 

P.idC.) = (25)! = "-'^ ^= (5) 

[ otherwise, 

where n,^r = 1]^^^,^,^. 

Note that apart from the boundary conditions, the weight W{C) is completely equivalent to that appears in the path 
integral representation of the system described by the Hamiltonian ( ^) . Therefore, for that part, we adopt the same 
labeling rule as the original S = ^ continuous- time loop algorithm [^|, which assigns a 'graph' Qw to a configuration 
C with labeling probability 

Tw{Qw\C) = . (6) 

Here, we follow the general framework of cluster algorithms presented in Ref. The weight Vw{Gw) is C-independent 
and non-negative, and AwiCyGw) is a compatibility function, which takes or 1. They satisfy 

W{C) ^ Vw{Gw)AwiC, Gw) (7) 

for any C. By this procedure the {d + l)-dimensional space is decomposed into a set of loops. Note that at this stage, 
some of the loops remain opened, because we have not defined any graphs for the boundary part. 

Next, for each local boundary weight Pi{dCi), we introduce (25*)! types of graphs {Gpi} {i = 1, ■ • ■ ,^s), each of 
which consists of 25* edges connecting one of subspins at r = to a subspin at r = /3 one by one (Fig. ||). We 
define a compatibility function Ap(9Ci, Gpi) as it takes 1 if every edge connects two subspins which have an identical 
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direction, or it takes otherwise. By using this compatibihty function, we can decompose the local boundary weight 
as 



P,{dC,) = gp,)/i2S)l , (8) 

because the number of compatible graphs for a configuration dCi is given by (25* — ni)lni\. We take the labeling 
probability for the local boundary weight as 

Apidc,,gp,)/{2sy. Ap{dc,,gp,) 

^^^^^^'^^^^ ^ mdC-) = (2S-n.)W ' 

that is, a graph is assigned out of the graphs compatible with dCi with equal probability. These boundary graphs make 
the remaining opened loops to be closed. By choosing the flipping probability to be free, that is, flipping spins on each 
loop simultaneously with probability ^ , one can show that the present stochastic process satisfies the detailed-balance 
condition 

The present algorithm includes the loop algorithm for the 5=1 antiferromagnetic Heisenberg model by Harada et 
aZ Q as a special case. However, as already seen, the present strategy does not depend so much on details of the model 
one considers; one can easily construct a cluster algorithm, if the mapped S — model has a cluster algorithm, which 
covers a model with general XYZ interaction or single-ion anisotropy JlQ], and the transverse- field Ising model [0. 



It can be applied also for a system with random size of spins and even for the classical Ising model [|2 13|. The 
resulting general-S" algorithm is ergodic, if the S — ^ cluster algorithm lying at the base is ergodic. The details 
of the algorithm for these models and the proof of its ergodicity will be presented elsewhere Note that the 

number of graphs introduced for the local boundary weight increases quite rapidly as S increases. However, the 
computational time in selecting a graph to be assigned is merely proportional to S, because the procedure is nothing 



but the random-permutation generation |14 



As an application of the algorithm, we simulated the antiferromagnetic Heisenberg chains with S" = 1, 2, and 3. It 
is conjectured by Haldane that the antiferromagnetic Heisenberg chain of integer spins has a finite excitation gap 
A(S') above its unique ground state, and the antiferromagnetic spin correlation along the chain decays exponentially 
with a finite correlation length £,x{S). For 5* = 1 and 2, a number of numerical studies have been accomplished 
(e.g., see [T^-pj|]) to confirm the validity of Haldane's conjecture. However, estimation of the first excitation gap of 
higher-spin chains has not yet been successful, since the magnitude of A(S') is considered to become exponentially 
small as S increases [ p^ . 

Consider a spin-S" chain of L sites at temperature T I//?)- We assume L is even. The correlation function of 
the staggered magnetization in the imaginary-time direction: 

Cir;L,P) = -^J2 fdt{{-lf-^\St{t)S^^{t + T)) (10) 

is an even function of t, and satisfies C(t -I- l3;L,j3) = C{t;L,P). At sufficiently low temperatures, the correlation 
function is expressed well as a sum of exponential functions, that is. 



C(r; L, (3) = ^ cosh 



T-p/2 



irAL) 



for /3 » Cr,o. (11) 



Here, we assume ^^- o > Cr.i > ^t,2 > • • • without loss of generality. The coefficients {ct} are directly related to the 
dynamic structure factor at momentum k = tt. In terms of ^.^,0(^)1 the gap to be estimated is given by 

A = lim — 5— . (12) 

In general, to solve Eq. ( pi] ) directly is extremely ill-posed [Q. However, as shown below, we can construct a 
systematic series of estimators at least for ^r,o(''")j if the coefficient Ci in Eq. ( pd| ) converges to zero rapidly enough 
for large i, and also if the difference between ^t- o(L) and ^T-^i(i) remains finite even in the thermodynamic limit. 
For S=l, these two conditions are numerically shown to be satisfied |^l[]. We expect similar situations for higher S, 
although there is no exact argument. 

For a given L, the well-known second- moment estimator of the correlation length, 
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FIG. 2. Temperature dependence of the spatial correlation length (open symbols) and the temporal correlation length 
(solid symbols) of the S = 1 (squares), 2 (triangles), and 3 (diamonds) antiferromagnetic Heisenberg chains. Statistical error 
of each data point is smaller than the line width. The system size is taken as L = 2S/T. 



converges to £^rfi{L) in the low-temperature limit, besides systematic corrections of 0{ai^r,i/^T,o) (* = 1, 2, • • ■). Here, 
C{uj) is the Fourier transform of the imaginary-time spin correlation function, i.e., (7(0;) = C{T)e^'^'^dT. In the 
present algorithm, C{uj) can be measured directly by means of the improved estimator |^^, which reduces the variance 
of the data greatly. We also consider a fourth-moment estimator; 




^~C^) = A. 33t^(^-l, (14) 



C(27r//3) - CiAn/P) 



which has smaller corrections of 0{ai{£^r.i/£,T.o)'^)- Construction of higher-order estimators is possible in a straight- 
forward way [|o| . 

The temperature dependence of S,-r{L,T) and also of the correlation length along the chain, S^x{L,T), is shown in 
Fig. ^. In the present simulation, the system size is taken as L = 2S/T. The integrated auto-correlation time of 
C(0) is of order unity and no significant sign of its growth is observed for larger L and 1/T. Measurement of physical 
quantities is performed after discarding first 10"^ Monte Carlo steps (MCS) for thermalization. One Monte Carlo step 
of the S = 3 chain with L = 5792 and T = 0.0010359, which is mapped to the system of 34752 subspins and 208512 
bonds prior to the simulation, takes about 16 seconds on 256 nodes of Hitachi SR-2201. 

As seen clearly in Fig ^, ^T-(L,r) and ^x{L,T) converge quite rapidly (probably exponentially) to a finite value 
for small T. We observe that for each S the data which satisfy the conditions, 1/T > 6^r and L > 6^x, exhibit no 
temperature (and system-size) dependence besides statistical fluctuations. The difference between the second-moment 
estimate (Eq. (|l3|)) and the fourth-moment one (Eq. (p^)) is invisible in the vertical scale of Fig ^. The difference 
between these two estimates at the lowest temperature is about 0.2% and 0.1% for S = 1 and S = 2, respectively. For 
S = 3, both coincide within the statistical errors. Furthermore, it is conflrmed that there is no significant difference 
between the fourth- moment estimate and the sixth- moment one even for S — 1. 



TABLE I. Ground-state energy density E/L, staggered susceptibility Xs, spatial correlation length (^x, and first excited gap 
A of the S = 1, 2, and 3 antiferromagnetic Heisenberg chains. Note that for each S, the physical quantities are estimated 
by a single Monte Carlo run on the system of size L at temperature T, which are presented in the second and third columns, 
respectively. 



s 


L 


T 


MCS 


E/L 






A 


1 


128 


0.0156250 


2 X 10' 


-1.401481(4) 


18.4048(7) 


6.0153(3) 


0.41048(6) 


2 


724 


0.0055249 


2 X 10^ 


-4.761249(6) 


1164.0(2) 


49.49(1) 


0.08917(4) 


3 


5792 


0.0010359 


3 X 10* 


-10.1239(1) 


158000(310) 


637(1) 


0.01002(3) 
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Thus, by using the fourth-moment estimator, we conclude 



( 0.41048(6) for 5* = 1 
A{S) = i 0.08917(4) for 5 = 2 (15) 
[ 0.01002(3) for 5 = 3 

as the magnitude of the Haldane gap. The results for other physical quantities, such as the energy density and the 
staggered susceptibility, are presented in Table ||. It should be emphasized that the present results are obtained 
without any extrapolation procedure; they are simply obtained by a single Monte Carlo run on the largest system at 
the lowest temperature for each S. 

For S — 1, the present estimate is completely consistent with 0.41050(2) and 0.41049(2) obtained by the DMRG 
calculation |l^ and by the exact diagonalization jl^], respectively. For 5 = 2, on the other hand, the numerical 
uncertainty in the present estimate is much smaller than in the previous studies (see, e.g., TABLE I in Ref. |p^ ). Fur- 
thermore, our estimate is slightly larger than 0.0876(13), which is obtained by the most recent DMRG calculation [p^ . 
In the DMRG study, for some technical reasons, open boundary conditions have been used, which is known to give 
quite large systematic corrections compared to the periodic boundary conditions. The reason of the disagreement 
might be due to an inappropriate scaling assumption in the DMRG study. Finally, as for the 5 = 3 case, it might be 
practically impossible to estimate the value of the gap by other numerical methods. The present result is completely 
new to our best knowledge. 
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